3D Genome Reconstruction from Partially Phased Hi-C Data

The 3-dimensional (3D) structure of the genome is of significant importance for many cellular processes. In this paper, we study the problem of reconstructing the 3D structure of chromosomes from Hi-C data of diploid organisms, which poses additional challenges compared to the better-studied haploid setting. With the help of techniques from algebraic geometry, we prove that a small amount of phased data is sufficient to ensure finite identifiability, both for noiseless and noisy data. In the light of these results, we propose a new 3D reconstruction method based on semidefinite programming, paired with numerical algebraic geometry and local optimization. The performance of this method is tested on several simulated datasets under different noise levels and with different amounts of phased data. We also apply it to a real dataset from mouse X chromosomes, and we are then able to recover previously known structural features.


Introduction
The eukaryotic chromatin has a three-dimensional (3D) structure in the cell nucleus, which has been shown to be important in regulating basic cellular functions, including gene regulation, transcription, replication, recombination, and DNA repair [43,45].The 3D DNA organization is also associated with brain development and function; in particular, it is shown to be misregulated in schizophrenia [33,35] and Alzheimer's disease [29].
All genetic material is stored in chromosomes, which interact in the cell nucleus, and the 3D chromatin structure influences the frequencies of such interactions.A benchmark tool to measure such frequencies is high-throughput chromosome conformation capture (Hi-C) [17].Hi-C first crosslinks cell genomes, which "freezes" contacts between DNA segments.Then the genome is cut into fragments, the fragments are ligated together and then are associated with equally-sized segments of the genome using high-throughput sequencing [34].These segments of the genome are called loci, and their size is known as resolution (e.g., bins of size 1Mb or 50Kb).The result of Hi-C is stored in a matrix called contact matrix whose elements are the contact counts between pairs of loci.
According to the structure they generate, computational methods for inferring the 3D chromatin structure from a contact matrix fall into two classes: ensemble and consensus methods.In a haploid setting (organisms having a single set of chromosomes), ensemble models such as MCMC5C [36], BACH-MIX [12] and Chrom3D [31], try to account for structure variations on the genome across cells by inferring a population of 3D structures.On the other hand, consensus methods aim at reconstructing one single 3D structure which may be used as a model for further analysis.In this category, probability-based methods such as PASTIS [44,5] and ASHIC [47] model contact counts as Poisson random variables of the Euclidean distances between loci, and distance-based methods such as ChromSDE [48] and ShRec3D [18] model contact counts as functions of the Euclidean distances.An extensive overview of different 3D genome reconstruction techniques is given in [30].
Most of the methods for 3D genome reconstructions from Hi-C data are for haploid organisms.However, like most mammals, humans are diploid organisms, in which the genetic information is stored in pairs of chromosomes called homologs.Homologous chromosomes are almost identical besides some single nucleotide polymorphisms (SNPs) [19].In the case of diploid organisms, the Hi-C data does not generally differentiate between homologous chromosomes.If we model each chromosome as a string of beads, then we associate two beads to each locus i ∈ {1, . . ., n}, one bead for each homolog.Therefore, each observed contact count c i,j between loci i and j represents aggregated contacts of four different types of interactions, more precisely one of the two homologous beads associated to locus i gets in contact with one of the two homologous beads associated to locus j, see Figure 1.This means that the Hi-C data is unphased.Phased Hi-C data that distinguishes contacts for homologs is rare.In our setting, we assume that the data is partially phased, i.e., some of the contact counts can be associated with a homolog.For example, in the (mouse) Patski (BL6xSpretus) [7,47] cell line, 35.6% of the contact counts are phased; while this value is as low as 0.14% in the human GM12878 cell line [34,47].Therefore, methods for inferring diploid 3D chromatin structure need to take into account the ambiguity of diploid Hi-C data to avoid inaccurate reconstructions.
Figure 1.Ambiguity of phased data.Each entry c i,j of the Hi-C matrix corresponds to four different contacts between the two pairs (x i , y i ) for locus i and (x j , y j ) for locus j.
Methods for 3D genome reconstruction in diploid organisms have been studied in [42,47,5,24,2,23,38].One approach is to phase Hi-C data [42,24,23], for example by assigning haplotypes to contacts based on assignments at neighboring contacts [42,23].Cauer et al. [5] and Ye and Ma [47] model contact counts as Poisson random variables.To find the optimal 3D chromatin structure, Cauer et al. maximize the associated likelihood function combined with two structural constraints.The first constraint imposes that the distances between neighboring beads are similar, and the second one requires that homologous chromosomes are located in different regions of the cell nucleus.On the other hand, Ye and Ma first compute the maximum likelihood estimate of model parameters for each of the homologs separately; these estimates are then refined by estimating the distance between the homologs.Belyaeva et al. [2] show identifiability of the 3D structure when the Euclidean distances between neighboring beads and higher-order contact counts between three or more loci simultaneously are given.Under these assumptions, the 3D reconstruction is obtained by combining distance geometry with semidefinite programming.Segal [38] applies recently developed imaging technology, in situ genome sequencing (IGS) [32], to point out issues in the assumptions made in [42,5,2], and suggests as alternative assumptions that intra-homolog distances are smaller than corresponding inter-homolog distances and intra-homolog distances are similar for homologous chromosomes.IGS [32] provides yet another method for inferring the 3D structure of the genome, however, at present the resolution and availability of IGS data is limited.

Contributions.
In this work, we focus on a distance-based approach for partially phased Hi-C data.In particular, we assume that contacts only for some loci are phased.In the string of beads model, the locations of the pair of beads associated to i-th loci are denoted by x i , y i ∈ R 3 .Then homologs are represented by two sequences x 1 , x 2 , . . ., x n and y 1 , x 2 , . . ., y n in R 3 ; see Figure 1.Inferring the 3D chromatin structure corresponds to estimating the bead coordinates.Based on Lieberman-Aiden et al. [22], we assume the power law dependency c i,j = γd α i,j , where α is a negative conversion factor, between the distance d i,j and contact count c i,j of loci i and j.Following Cauer et al. [5], we assume that a contact count between loci is given by the sum of all possible contact counts between the corresponding beads.We call a bead unambiguous if the contacts for the corresponding locus are phased; otherwise, we call a bead ambiguous.
Our first main contribution is to show that for negative rational conversion factors α, knowing the locations of six unambiguous beads ensures that there are generically finitely many possible locations for the other beads, both in the noiseless (Theorem 3.1) and noisy (Corollary 3.5) setting.Moreover, we prove finite identifiability also in the fully ambiguous setting when α = −2 and the number of loci is at least 12 (Theorem 3.6).Note that the identifiability does not hold for α = 2 as shown in [2].
Our second main contribution is to provide a reconstruction method when α = −2, based on semidefinite programming combined with numerical algebraic geometry and local optimization (section 4).The general idea is the following: We first estimate the coordinates of the unambiguous beads using only the unambiguous contact counts (which precisely corresponds to the haploid setting) using the SDP-based solver implemented in ChromSDE [48].We then exploit our theoretical result on finite identifiability to estimate the coordinates of the ambiguous beads, one by one, by solving several polynomial systems numerically.These estimates are then improved by a local estimation step considering all contact counts.Finally, a clustering algorithm is used to overcome the symmetry (x i , y i ) → (y i , x i ) in the estimation for the ambiguous beads.
The paper is organized as follows.In section 2, we introduce our mathematical model for the 3D genome reconstruction problem.In section 3, we recall identifiability results in the unambiguous setting (section 3.1) and then prove identifiability results in the partially ambiguous setting (section 3.2) and in the fully ambiguous setting (section 3.3).We describe our reconstruction method in section 4. We test the performance of our method on synthetic datasets and on a real dataset from the mouse X chromosomes in section 5. We conclude with a discussion about future research directions in section 6.

Mathematical model for 3D genome reconstruction
In this section we introduce the distance-based model under which we study 3D genome reconstruction.In section 2.1 we give the background on contact count matrices.In section 2.2 we describe a power-law between contacts and distances, which allows to translate the information about contacts into distances.
2.1.Contact count matrices.We model the genome as a string of 2n beads, corresponding to n pairs of homologous beads.The positions of the beads are recorded by a matrix The positions x i and y i correspond to homologous beads.When convenient, we use the notation z 1 := x 1 , . . ., z n := x n , z n+1 := y 1 , . . ., z 2n := y n .In this notation, Let U be the subset of pairs that are unambiguous, i.e., beads in the pair can be distinguished, and let A be the subset of pairs that are ambiguous, i.e., beads in the pair cannot be distinguished.The sets U and A form a partition of [n].
A Hi-C matrix C is a matrix with each row and column corresponding to a genomic locus.Following Cauer et al. [5], we call these contact counts ambiguous and denote the corresponding contact count matrix by C A .If parental genotypes are available, then one can use SNPs to map some reads to each haplotype [7,25,34].If both ends of a read contains SNPs that can be associated to a single parent, then the contact count is called unambiguous and the corresponding contact count matrix is denoted by C U .Finally, if only one of the genomic loci present in an interaction can be mapped to one of the homologous chromosomes, then the count is called partially ambiguous and the contact count matrix is denoted by C P .
The unambiguous count matrix C U is a 2n × 2n matrix with the first n indices corresponding to x 1 , . . ., x n and the last n indices corresponding to y 1 , . . ., y n .The ambiguous count matrix C A is an n×n matrix and we assume that each ambiguous count is the sum of four unambiguous counts: The partially ambiguous count matrix C P is a 2n×n matrix and each partially ambiguous count is the sum of two unambiguous counts: Seven different types of contacts between the ith and jth locus.

Contacts and distances.
Denoting the distance ∥z i − z j ∥ between z i and z j by d i,j , the power law dependency observed by Lieberman-Aiden et al. [22] can be written as where α < 0 is a conversion factor and γ > 0 is a scaling factor.This relationship between contact counts and distances is assumed in [2,48], while in [5,44] the contact counts c i,j are modeled as Poisson random variables with the Poisson parameter being βd α i,j .
In our paper, we assume that contact counts are related to distances by (2.1).Similarly to [2], we set γ = 1 and in parts of the article α = −2.In general, the conversion factor α depends on a dataset and its estimation can be part of the reconstruction problem [44,48].Setting γ = 1 means that we recover the configuration up to a scaling factor.In practice, the configuration can be rescaled using biological knowledge, e.g., the radius of the nucleus.
Our approach to 3D genome reconstruction builds on the power law dependency between contacts and distances between unambiguous beads.We convert the empirical contact counts to Euclidean distances and then aim to reconstruct the positions of beads from the distances.This leads us to the following system of equations: If α is an even integer, then (2.2) is a system of rational equations.
Determining the points x i , y i , where i ∈ U , is the classical Euclidean distance problem: We know the (noisy) pairwise distances between points and would like to construct the locations of points, see section 3.1 for details.Hence after section 3.1 we assume that we have estimated the locations of points x i , y i , where i ∈ U , and we would like to determine the points x i , y i , where i ∈ A.

Identifiability
In this section, we study the uniqueness of the solutions of the system (2.2) up to rigid transformations (translations, rotations and reflections), or in other words, the identifiability of the locations of beads.We study the unambiguous, partially ambiguous and ambiguous settings in sections 3.1, 3.2 and 3.3, respectively.
3.1.Unambiguous setting and Euclidean distance geometry.If all pairs are unambiguous, i.e., U = [n], then constructing the original points translates to a classical problem in Euclidean distance geometry.The principal task in Euclidean distance geometry is to construct original points from pairwise distances between them.In the rest of the subsection, we will recall how to solve this problem.Since pairwise distances are invariant under translations, rotations and reflections (rigid transformations), then the original points can be reconstructed up to rigid transformations.For an overview of distance geometry and Euclidean distance matrices, we refer the reader to [8,16,21,27].
The Gram matrix of the points z 1 , . . ., z 2n is defined as i=1 z i and zi = z i − z for i = 1, . . ., 2n.The matrix Z = [z 1 , . . ., z2n ] T gives the locations of points after centering them around the origin.Let G denote the Gram matrix of the centered point configuration z1 , . . ., z2n .
Let D i,j = ∥z i −z j ∥ 2 denote the squared Euclidean distance between the points z i and z j .The Euclidean distance matrix of the points z 1 , . . ., z 2n is defined as D = (D i,j ) 1≤i,j≤2n ∈ R 2n×2n .To express the centered Gram matrix in terms of the Euclidean distance matrix, we define the geometric centering matrix where I 2n is the 2n × 2n identity matrix and 1 is the vector of ones.The linear relationship between G and D is given by G = − 1 2 JDJ.
Therefore, given the Euclidean distance matrix, we can construct the centered Gram matrix for the points z 1 , . . ., z 2n .
The centered points up to rigid transformations are extracted from the centered Gram matrix G using the eigendecomposition G = QΛQ −1 , where Q is orthonormal and Λ is a diagonal matrix with entries ordered in decreasing order 3 .In the case of noiseless distance matrix D, the Gram matrix G has rank three and the diagonal matrix Λ has precisely three non-zero entries.Hence we could obtain Ẑ also from QΛ 1/2 by truncating zero columns.Using Λ has the advantage that it gives an approximation for the points also for a noisy distance matrix D. The uniqueness of Ẑ up to rotations and reflections follows from [15, Proposition 3.2], which states that AA T = BB T if and only if A = BQ for some orthogonal matrix Q.
The procedure that transforms the distance matrix to origin centered Gram matrix and then uses eigendecomposition for constructing original points is called classical multidimensional scaling (cMDS) [6].Although cMDS is widely used in practice, it does not always find the distance matrix that minimizes the Frobenius norm to the empirical noisy distance matrix [40].Other approaches to solving the Euclidean distance and Euclidean completion problems include nonconvex [10,26] as well semidefinite formulations [1,11,28,46,48,49].

3.2.
Partially ambiguous setting.The next theorem establishes the uniqueness of the solutions of the system (2.2) in the presence of ambiguous pairs.In particular, it states that there are finitely many possible locations for beads in one ambiguous pair given the locations of six unambiguous beads.The identifiability results in this subsection hold for all negative rational numbers α.In the rest of the paper, we denote the true but unknown coordinates by x * and the symbol x stands for a variable that we want to solve for.We write ∥ • ∥ for the standard inner product on R 3 .Theorem 3.1.Let α be a negative rational number.Then for a * , b * , . . ., f * , x * , y * ∈ R 3 sufficiently general, the system of six equations in the six unknowns x 1 , x 2 , x 3 , y 1 , y 2 , y 3 ∈ R has only finitely many solutions.
Remark 3.2.The proof will show that this system has only finitely many solutions over the complex numbers.
We believe that the theorem holds for general nonzero rational α.Indeed, our argument works, with a minor modification, also for α > 2, but for α in the range (0, 2] a refinement of the argument is needed. Write α 2 = m n with m, n relatively prime integers, m ̸ = 0, and n > 0. Consider the affine variety X ⊆ (C 3 ) 8 × (C 2 ) 6 consisting of all tuples Note that, if x * , t * are real, then it follows that and similarly for Q(y * − t * ).Hence if a * , . . ., y * are all real, then the point is a point in X with real-valued coordinates.
The projection π from X to the open affine subset U ⊆ (C 3 ) 8 where all Q(x * −t * ) and Q(y * −t * ) are nonzero is a finite morphism with fibers of cardinality n 12 ; to see this cardinality note that there are n possible choices for each of the numbers r t * , s t * .Each irreducible component of X is a smooth variety of dimension 24.
Consider the map ψ : X → (C 3 × C 1 ) 6 defined by We claim that for q in some open dense subset of X, the derivative d q ψ has full rank 24.For this, it suffices to find one point p ∈ U such that d q ψ has rank 24 at each of the n 12 points q ∈ π −1 (p).We take a real-valued point p := (a * , b * , . . ., f * , x * , y * ) ∈ (R 3 ) 8 to be specified later on.Let q ∈ π −1 (p).Then, near q, the map ψ factorises via π and the unique algebraic map where ξ t * and ζ t * are n-th roots of unity in C depending on which q is chosen among the n 12 points in π −1 (p).The situation is summarised in the following diagram: and since d q π is a linear isomorphism, it suffices to prove that d p ψ ′ is a linear isomorphism.Suppose that (a ′ , . . ., f ′ , x ′ , y ′ ) ∈ ker d p ψ ′ .Then, since the map ψ ′ remembers a, . . ., f , it follows immediately that a ′ = . . .= f ′ = 0. On the other hand, by differentiating we find that, for each t * ∈ {a * , . . ., f * }, where ⟨•, •⟩ stands for the standard bilinear form on C 3 .In other words, the vector (x ′ , y ′ ) ∈ C 6 is in the kernel of the 6 × 6-matrix where we have interpreted a * , . . ., f * , x * , y * as row vectors.It suffices to show that, for some specific choice of p = (a * , . . ., f * , x * , y * ) ∈ (R 3 ) 8 , this matrix is nonsingular for all n 12 choices of ((ξ t * , η t * )) t * .
Then the matrix M becomes, with β = α − 2: where R is a sum of (products of) roots of unity.Now α < 0 implies that β < −2, so that 2 + 3β < 2β < 0. Since roots of unity have 2-adic valuation 0, the second term in the expression above is the unique term with minimal 2-adic valuation.Hence det(M ) ̸ = 0, as desired.
It follows that ψ is a dominant morphism from each irreducible component of X into (C 3 × C 1 ) 6 , and hence for all q in an open dense subset of X, the fiber ψ −1 (ψ(q)) is finite.This then holds, in particular, for q in an open dense subset of the real points as in (3.2).This proves the theorem.□ Remark 3.3.If α > 2, then β > 0, and hence the unique term with minimal 2-adic valuation in (3.3) is the first term.This can be used to show that the theorem holds then, as well.The only subtlety is that for positive α, solutions where x or y equal one of the points a * , . . ., f * are not automatically excluded, and these are not seen by the variety X.But a straightforward argument shows that such solutions do not exist for sufficiently general choices of a * , . . ., f * , x * , y * .
We now consider the setting when we know locations of seven unambiguous beads.In the special case when α = −2, we construct the ideal generated by the polynomials obtained from rational equations (3.1) for seven unambiguous beads after moving all terms to one side and clearing the denominators.Based on symbolic computations in Macaulay2 for the degree of this ideal, we conjecture that the location of a seventh unambiguous bead guarantees unique identifiability of an ambiguous pair of beads: Conjecture 3.4.Let a * , b * , c * , d * , e * , f * , g * , x * , y * ∈ R 3 be sufficiently general.The system of rational equations has precisely two solutions (x * , y * ) and (y * , x * ).
In practice, we only have noisy estimates a, b, . . ., f ∈ R 3 of the true positions of unambiguous beads a * , b * , . . ., f * ∈ R 3 , and we have noisy observations c t of the true contact counts c * t := ∥x * − t * ∥ α + ∥y * − t * ∥ α .We aim to find x, y ∈ R 3 such that ∥x − t∥ α + ∥y − t∥ α = c t for t = a, b, . . ., f.We may write c t = ∥x * − t∥ α + ∥y * − t∥ α + ϵ t for some ϵ t that depends on the noise level.Hence, the above system of equations can be rephrased as ∥x − t∥ α + ∥y − t∥ α = ∥x * − t∥ α + ∥y * − t∥ α + ϵ t for t = a, b, . . ., f. (3.5) In the following corollary we show that this system has generically finitely many solutions.
We showed that ψ is a dominant morphism from each irreducible component of X into (C 3 ×C 1 ) 6 , and that each irreducible component of X is 24-dimensional.Every solution to (3.6) is the (x, y)component of a point in the fiber Since this is a fiber over a sufficiently general point, the fiber is finite.□ Corollary 3.5 will be the basis of a numerical algebraic geometric based reconstruction method in section 4.
3.3.Ambiguous setting.Finally we consider the ambiguous setting, where one would like to reconstruct the locations of beads only from ambiguous contact counts.It is shown in [2] that for α = 2, one does not have finite identifiability no matter how many pairs of ambiguous beads one considers.We show finite identifiability for the locations of beads given contact counts for 12 pairs of ambiguous beads for α = −2 in both the noisy and noiseless setting.We believe that the result might be true for further conversion factors α's, however our proof technique does not directly generalize.Theorem 3.6.Let α = −2.Then for (c ij ) 1≤i<j≤12 ∈ R 66 sufficiently general, the system of 66 equations in the 72 unknowns x 1,1 , x 1,2 , x 1,3 , y ) ∈ (R 3 ) 24 , the system has finitely many solutions up to rigid transformation.
Proof.As before, we write Q(x . By a computer calculation (with exact arithmetic) we found that at a randomly chosen q ∈ X with rational coordinates, the derivative d q ψ had full rank 66.It then follows that for q in some open dense subset of X, d q ψ has rank 66. Hence ψ is dominant, and for any sufficiently general c ∈ C 66 , all irreducible components of the fiber ψ −1 (c) have dimension 6.Moreover, each such component C is preserved by the 6-dimensional connected group The stabilizer in G of a sufficiently general point in X is zero-dimensional.This follows from a Lie algebra argument: if a point (x * 1 , y * 1 , . . ., x * 12 , y * 12 ) ∈ X has a positive-dimensional stabilizer in G, then there is a nonzero element A in the Lie algebra of SO(3, C) that maps all the differences x * i − x * j , x * i − y * j , y * i − y * j to zero.Since A is a skew-symmetric matrix and hence of rank 2, it follows that all points x * i , y * j lie on a line.The variety of such collinear tuples has dimension 28, so it does not map dominantly to C 66 .Hence there exists a Zariski open dense subset V ⊆ C 66 such that for all c ∈ V , the fiber ψ −1 (c) contains no points with positive-dimensional stabilizers in G, and hence ψ −1 (c) is a disjoint union of finitely many 6-dimensional G-orbits.Likewise, ψ −1 (V ) is a Zariski open dense subset of (C 3 ) 24 such that ψ −1 (ψ(q)) consists of finitely many G-orbits for all q ∈ ψ −1 (V ).With this, we have proven the complex analog of the theorem.
To obtain the statement over the real numbers, we note that if c ∈ V has real-valued coordinates, then a finite number of the G-orbits that make up ψ −1 (c) contain a real-valued tuple.If G • q for q ∈ (R 3 ) 24 is such an orbit, it holds that (G • q) ∩ (R 3 ) 24 = (SO(3, R) ⋉ R 3 ) • q whenever the 24 points that make up the tuple q are not coplanar.The set of coplanar configurations form a subset of X of dimension 51, and does therefore not map dominantly to C 66 .Hence, by shrinking V appropriately, we can assume that no fibers above it contain coplanar configurations.In particular, this means that the real part of the fiber over any real point in V consists of a finitely many orbits under the action of SO(3, R) ⋉ R 3 , as desired.□ Remark 3.7.A standard numerical algebraic geometry computation with monodromy and the certification techniques of [3], using HomotopyContinuation.jl (see, e.g., [41]), proves that the system (3.6)generically has more than 1000 complex solutions up to the action of O(3, C) ⋉ C 3 and the symmetries (x i , y i ) → (y i , x i ) for i = 1, . . ., 12.This constitutes theoretical motivation for working with partially phased data, even if we, in principle, have finite identifiability already from the unphased data.
Remark 3.8.When α = 2, which corresponds to the setting studied in [2], then computationally we found that for some special choices of x * 1 , y * 1 , . . ., x * 12 , y * 12 ∈ R 3 the rank of the Jacobian matrix in Theorem 3.6 is 42.This is consistent with the fact that Theorem 3.6 fails for α = 2 [2].

A new reconstruction method
In this section, we outline a new approach to diploid 3D genome reconstruction for partially phased data, based on the theoretical results discussed in subsection 3.2.The method consists of the following main steps: (1) Estimation of the unambiguous beads {x i , y i } i∈U through semidefinite programming (discussed in subsection 4.1).(2) A preliminary estimation of the ambiguous beads using numerical algebraic geometry, based on Corollary 3.5 (discussed in subsection 4.2).(3) A refinement of this estimation using local optimization (discussed in subsection 4.3).(4) A final clustering step, where we disambiguate between the estimations (x i , y i ) and (y i , x i ) for each i ∈ A, based on the assumption that homolog chromosomes are separated in space (discussed in subsection 4.4).
In what follows, we will refer to this method by the acronym SNLC (formed from the initial letters in semidefinite programming, numerical algebraic geometry, local optimization and clustering).
4.1.Estimation of the positions of unambiguous beads.As discussed in section 3.1, the unambiguous bead coordinates {x i , y i } i∈U = {z i } i∈U ∪(n+U ) can be estimated with semidefinite programming.More specifically, we use ChromSDE [48, Section 2.1] for this part of our reconstruction, which relies on a specialized solver from [14], to solve an SDP relaxation of the optimization problem min with λ = 0.01 (cf.[48,Equation 4]).The terms in the first sum are weighted by the square root for the corresponding contact counts, in order to account for the fact that higher counts can be assumed to be less susceptible to noise.
4.2.Preliminary estimation using numerical algebraic geometry.To estimate the coordinates of the ambiguous beads {x i , y i } i∈A , we will use a method based on numerical equation solving, where we estimate the ambiguous bead pairs one by one.
Let x, y be the unknown coordinates in R 3 of a pair of ambiguous beads.We pick six unambiguous beads with already estimated coordinates a, b, c, d, e, f ∈ R 3 .For each t ∈ {a, . . ., f }, let c t ∈ R be the corresponding partially ambiguous counts between t and the ambiguous bead pair (x, y).Clearing the denominators in the system (3.6), we obtain a system of polynomial equations ∥x − t∥ 2 + ∥y − t∥ 2 = c t ∥x − t∥ 2 ∥y − t∥ 2 for t = a, b, c, d, e, f. ( By Corollary 3.5, this system has finitely many complex solutions both in the noiseless and noisy setting, which can be found using homotopy continuation. We observe that the system (4.2) generally has 80 complex solutions, and we only expect one pair of solutions (x, y), (y, x) to correspond to an accurate estimation.Naively adding another polynomial arising from a seventh unambiguous bead (as in Conjecture 3.4) does not work; in the noisy setting this over-determined system typically lacks solutions.Instead, we compute an estimation based on the following two heuristic assumptions: (1) The most accurate estimation should be approximately real, in the sense that the maxnorm of the imaginary part is below a certain tolerance (in this work, 0.15 was used for the experiments in both subsections 5.1 and 5.2).The choice of this threshold was made based on analysing the imaginary parts of solutions to (4.2) for various choices of unambiguous beads, see Figure 9. (2) The most accurate estimation should be consistent when we change the choice of six unambiguous beads.
Based on these assumptions, we apply the following strategy.We make a number N ≥ 2, choices of sets of six unambiguous beads, and solve the corresponding N square systems of the form (4.2).Since larger contact counts can be expected to have smaller relative noise, we make the choices of beads among the 20 unambiguous beads t that have highest contact count c t to the ambiguous locus at hand.For each system, we pick out the approximately real solutions, and obtain N sets S 1 , . . ., S N ⊆ R 6 consisting of the real parts of the approximately real solutions.Up to the symmetry (x, y) → (y, x), we expect these sets to have a unique "approximately common" element.We therefore compute, by an exhaustive search, the tuple (w 1 , . . ., w N ) ∈ S 1 × • • • × S N that minimizes the sum and use as our estimation of (x, y).For the computations presented in section 5, we use N = 5.
To solve the systems, we use the Julia package HomotopyContinuation.jl [4], and follow the two-phase procedure described in [39,Section 7.2].For the first phase, we solve (4.2) with randomly chosen parameters a * , . . ., f * ∈ C 3 and c a * , . . ., c f * ∈ C, using a polyhedral start system [13].We trace 1280 paths in this first phase, since the Newton polytopes of the polynomials appearing in the system (4.2) all contain the origin, and have a mixed volume of 1280, which makes 1280 an upper bound on the number of complex solutions by [20,Theorem 2.4].For the second phase, we use a straight-line homotopy in parameter space from the randomly chosen parameters a * , . . ., f * ∈ C 3 and c a * , . . ., c f * ∈ C, to the values a, . . ., f and c a , . . ., c f ∈ C at hand.We observe that we generally find 80 complex solutions in the first phase, which means 40 orbits with respect to the symmetry (x, y) → (y, x).By the discussion in [39,Section 7.6], it is enough to only trace one path per orbit, so in the end, we only trace 40 paths in the second phase.
Remark 4.1.If the noise levels are sufficiently high, there could be choices of six unambiguous beads for which the system lacks approximately-real solutions.If this situation is encountered, we try to redraw the six unambiguous beads until we find an approximately-real solution.If this does not succeed within a certain number of attempts (100 in the experiments conducted for this paper), we use the average of the closest neighboring unambiguous beads instead.

Local optimization.
A disadvantage of the numerical algebraic geometry based estimation discussed in the previous subsection is that it only takes into account "local" information about the interactions for one ambiguous locus at a time, which might make it more sensitive to noise.In our proposed method, we therefore refine this preliminary estimation of {x i , y i } i∈A further in a local optimization step that takes into account the "global" information of all available data.
The idea is to estimate {x i , y i } i∈A by solving the optimization problem while keeping the estimates of {x i , y i } i∈U from the ChromSDE step fixed.We use the quasi-Newton method for unconstrained optimization implemented in the Matlab Optimization Toolbox for this step.The already estimated coordinates of {x i , y i } i∈A from the numerical algebraic geometry step are used for the initialization.

4.4.
Clustering to break symmetry.Our objective function remains invariant if we exchange x i and y i for any i ∈ A. We can break symmetry by relying on the empirical observation that homologous chromosomes typically are spatially separated in different so-called compartments of the nucleus [9].Let (x i , ȳi ) n i=1 denote the estimates from the previous steps.Our final estimations will be obtained by solving the minimization problem where (x i , y i ) = (x i , ȳi ) for i ∈ U are fixed, and (x i , y i ) ∈ {(x i , ȳi ), (ȳ i , xi )} for i ∈ A are the optimization variables.The optimal solution can be computed efficiently, as explained next.
We first decompose the problem into contiguous chunks of ambiguous beads.Let (i 1 , . . ., i L ) := U be the indices of the unambiguous beads and let i 0 := 1, i L+1 := n.The optimization problem can be phrased as where there is one summand G ℓ (x, y) for each contiguous chunk of ambiguous beads.Since the summands G ℓ (x, y) do not share any ambiguous bead, we can minimize them independently.
We proceed to describe the optimal solution of the problem.Let The variable s i indicates whether we keep using (x i , ȳi ) or we reverse it.Note that s i = 1 for i ∈ U .The next lemma gives the optimal assignment of s i for i ∈ A. This assignment is constructed by using inner products w i,i+1 .
Lemma 4.2.The optimal solution of (4.4) can be constructed as follows: (1) For the last chunk (ℓ = L) we have where sgn(•) is the sign function and sgn(0) can be either 1 or −1.
(2) For the first chunk (ℓ = 0) we have (3) For any other chunk, let k be the index of the smallest absolute value Proof.Denoting ūi := 1 2 (x i + ȳi ), vi := 1 2 (x i − ȳi ), then Since xi , ȳi , ūi , vi are constants, minimizing g i,i+1 (x, y) is equivalent to maximizing w i,i+1 s i s i+1 .Then for each chunk we have to solve the optimization problem max The formulas from the first and last chunk are such that w i,i+1 s * i s * i+1 ≥ 0 for all i.This is possible because in these cases only one of the endpoints has a fixed value, and the remaining values are computed recursively starting from such a fixed point.Since all summands are nonnegative, the sum in (4.6) is maximized.
For the inner chunks, the two endpoints are fixed, so it may not be possible to have that w i,i+1 s * i s * i+1 ≥ 0 for all indices.In an optimal assignment we should pick at most one term to be negative, and such a term (if it exists) should be the one with the smallest absolute value |w i,i+1 |.This leads to the formula from the lemma.□

Experiments
In this section, we apply the SNLC scheme described in section 4 to synthetic and real datasets, and compare its performance with the preexisting software packages ASHIC [47] and PASTIS [5].We chose these two reconstruction methods for comparison because they are best suited for our setting.Also Belyaeva et al. [2] and Tan et al. [42] have reconstruction methods for diploid organisms, but the former method requires higher-order contact information and the latter method is targeted for single cell data.
All SNLC experiments are done using Julia 1.6.1, with ChromSDE being run in Matlab 2021a, and the Julia package MATLAB.jl(v0.8.3) acting as interface between Julia and Matlab.The numerical algebraic geometry part of the estimation procedure is done with HomotopyContinuation.jl (v2.5.5) [4].The PASTIS experiments are run in Python 3.8.10,and the ASHIC experiments in Python 3.10.5.
For the PASTIS computations, we fix α = −2 to ensure compatibility with the modelling assumptions made in this paper.We run PASTIS without filtering, in order to make it possible to compare RMSD values.Since PASTIS only takes integer inputs, we multiply the theoretical contact counts calculated by (2.2) by a factor 10 5 and round them to the nearest integer.Following the approach taken in [5], we use a coarse grid search to find the optimal coefficients for the homolog separating constraint and bead connectivity constraints.Specifically, we fix a structure simulated with the same method as used in the experiments, and compute the RMSD values for all λ 1 , λ 2 ∈ {1, 10 1 , 10 2 , . . ., 10 12 }.In this way, we find that λ 1 = 10 11 and λ 2 = 10 12 give optimal results.
For the ASHIC computations, we use the ASHIC-ZIPM method, which has the lowest distance error rate among the ASHIC's models according to [47, Figure 2] and models the contact counts as a zero-inflated Poisson distribution (ZIP) to account for the sparsity of the Hi-C matrix.We run ASHIC without filtering out any loci and with the setting aggregate to ensure that the coordinates of all beads are estimated.
5.1.Synthetic data.We conduct a number of experiments where we simulate a single chromosome pair (referred to as X and Y in figures) through Brownian motion with fixed step length, compute unambiguous, partially ambiguous and ambiguous contact counts according to (2.2), add noise, and then try to recover the structure of the chromosomes through the SNLC scheme described in section 4. Following [2], we model noise by multiplying each entry of C U , C P and C A by a factor 1 + δ, where δ is sampled uniformly from the interval (−ε, ε) for some chosen noise level ε ∈ [0, 1].
As a measure of the quality of the reconstruction, we use the minimal root-mean square distance (RMSD) between, on the one hand, the true coordinates (x * i , y * i ) n i=1 , and, on the other hand, the estimated coordinates (x i , y i ) n i=1 after rigid transformations and scaling, i.e., we find the minimum min This can be seen as a version of the classical Procrustes problem solved in [37], which is implemented in Matlab as the function procrustes.
Specific examples of reconstructions of the Brownian motion and helix-shaped chromosomes obtained with SNLC at varying noise levels and 50% of ambiguous beads are shown in Figure 3.For low noise levels the reconstructions by SNLC and the original structure highly overlap.For higher noise levels the general region occupied by the reconstructions overlaps with the original structure, while the local features become less aligned.Analogous reconstructions obtained with SNLC without the local optimization step are shown in Figure 6 in Appendix.
A comparison of how the quality of the reconstruction depends on the noise level and proportion of ambiguous beads for SNLC, ASHIC and PASTIS is done in Figure 4. We measure the RMSD value between the reconstructed and original 3D structure for different noise levels over 20 runs.The RMSD values obtained by SNLC are consistently lower than the ones obtained by ASHIC and PASTIS.The difference is specially large for low to medium noise levels.While our method outperforms ASHIC and PASTIS in the setting considered in this paper, it is worth mentioning that ASHIC and PASTIS work also in a more general setting, where there might be contacts of all three types (ambiguous, partially ambiguous and unambiguous) between every pair of loci.5.2.Experimentally obtained data.We compute SNLC reconstructions based on the real dataset explored in [5], which is obtained from Hi-C experiments on the X chromosomes in the Patski (BL6xSpretus) cell line.The data has been recorded at a resolution of 500 kb, which corresponds to 343 bead pairs in our model.
For some of these pairs, no or only very low contact counts have been recorded.Since such low contact counts are susceptible to high uncertainty and can be assumed to be a consequence of experimental errors, we exclude the 47 loci with the lowest total contact counts from the analysis.To select the cutoff, the loci are sorted according to the total contact counts (see Figure 7 (a) in Appendix), and the ratios between the total contact counts for consecutive loci are computed.A peak for these ratios is observed at the 47th contact count, as shown in Figure 7 (b) in Appendix.After applying this filter, we obtain a dataset with 296 loci.Out of these, we consider as ambiguous all loci i for which less than 40% of the total contact count comes from contacts where x i and y i were not distinguishable.These proportions for all loci are shown in Figure 7 (c) in Appendix.For the Patski dataset, we obtain 46 ambiguous loci and 250 unambiguous loci in this way.In the Patski dataset, a locus can simultaneously participate in unambiguous, partially ambiguous and ambiguous contacts.To obtain the setting of our paper where loci are partitioned into unambiguous or ambiguous, we reassign the contacts according to whether a locus is unambiguous or ambiguous.Our reassignment method is motivated by the assignment of haplotype to unphased Hi-C reads in [23].The exact formulas are given in Appendix.
The reconstruction obtained via SNLC can be found in Figure 5 (a).The logarithmic heatmaps for contact count matrices for original data and the SNLC reconstruction are shown in Figure 8.
It was discovered in [7] that the inactive homolog in the Patski X chromosome pair has a bipartite structure, consisting of two superdomains with frequent intra-chromosome contacts within the superdomains and a boundary region between the two superdomains.The active homolog does not exhibit the same behaviour.The boundary region on the inactive X chromosome is centered at 72.8-72.9MB [7] which at the 500 kB resolution corresponds to the bead 146 [5].We show in Figure 5 (b) that the two chromosomes reconstructed using SNLC exhibit this structure by computing the bipartite index for the respective homologs as in [5,7].We recall that, in the setting of a single chromosome with beads z 1 , . . ., z n ∈ R 3 , the bipartite index is defined as the ratio of intra-superdomain to inter-superdomain contacts in the reconstruction:

Discussion
In this article we study the finite identifiability of 3D genome reconstruction from contact counts under the model where the distances d i,j and contact counts c i,j between two beads i and j follow the power law dependency c i,j = d α i,j for a conversion factor α < 0. We show that if at least six beads are unambiguous, then the locations of the rest of the beads can be finitely identified from partially ambiguous contact counts for rational α satisfying α < 0 or α > 2.
In the fully ambiguous setting, we prove finite identifiability for α = −2, given ambiguous contact counts for at least 12 pairs of beads.From [2] it is known that finite identifiability does not hold in the fully ambiguous setting for α = 2.It is an open question whether finite identifiability of 3D genome reconstruction holds for other α ∈ R\{−2, 2} in the fully ambiguous setting and for rational α ∈ (0, 2] in the partially ambiguous setting.We conjecture that in the partially ambiguous setting seven unambiguous loci guarantee unique identifiability of the 3D reconstruction for rational α < 0 or α > 2. When α = −2, then one approach to studying the unique identifiability might be via the degree of a parametrized family of algebraic varieties.
After establishing the identifiability, we suggest a reconstruction method for the partially ambiguous setting with α = −2 that combines semidefinite programming, homotopy continuation in numerical algebraic geometry, local optimization and clustering.To speed up the homotopy continuation based part, we observe that the parametrized system of polynomial equations corresponding to six unambiguous beads has 40 pairs of complex solutions and we trace one path for each orbit.It is an open question to prove that for sufficiently general parameters the system has 40 pairs of complex solution.This question again reduces to studying the degree of a family of algebraic varieties.While our goal is to highlight the potential of our method, one could further regularize its output and use interpolation for the beads that are far away from the neighboring beads.A future research direction is to explore whether numerical algebraic geometry or semidefinite programming based methods can be proposed also for other conversion factors α < 0.

Appendix
In this part of the paper, we include additional details and figures for the experiments in section 5.
Figure 6 shows reconstructions of the same chromosomes as displayed in Figure 3 but without the local optimization step, indicating that semidefinite programming, numerical algebraic geometry and clustering alone can recover the main features of the 3D structure.Figure 7 illustrates the preprocessing steps of the real dataset where loci with low contact counts are removed and the rest of the loci are partitioned into unambiguous and ambiguous.The total contact count for the ith locus is defined as the sum of all contacts where it participates: c A (i, j)+ c P (i, j)+ c P (i + n, j) + j∈[2n] c P (j, i)+c U (i, j)+c U (i + n, j) .
Similarly, we define the unambiguity quotient as the proportion of T (i) that consists of contacts where x i and y i could be distinguished: The i-th smallest total contact count The peak corresponding to the ratio between the 48th and the 47th smallest count is used as a motivation for excluding the 47 loci with smallest total contact from the analysis.(c) Unambiguity quotients for each of the remaining 296 loci, sorted in increasing order.We consider a locus as ambiguous if this ratio is less than 0.4; otherwise, we consider it as unambiguous.
To obtain the setting of our paper where loci are partitioned into unambiguous or ambiguous, we reassign the contact counts of CU CP and CA of the Patski dataset according to whether a locus unambiguous or ambiguous.For i, j ∈ U , we define c U i,j = cU i,j + cP .
Finally, for i, j ∈ A, we define c A i,j = cU i,j + cU i,j+n + cU i+n,j + cU i+n,j+n + cP i,j + cP i+n,j + cP j,i + cP j+n,i + cA i,j .
In Figure 8 in Appendix, the experimental contact counts from the Patski dataset are compared with the contact counts from the SNLC reconstruction.
Figure 9 shows how the max-norm of the imaginary part of the solutions varies between different instances of the system (4.2) used for the reconstruction in Figure 3(b), and for the reconstruction from the Patski data in Figure 5.A complete set of figures for these two datasets can be found in the Github repository.Taken together, the figures indicate that a max-norm of 0.15 was an appropriate threshold for approximate realness for both data sets, in the sense that it is low enough to single out solutions that have significantly smaller imaginary parts than the others, while also ensuring that it is possible to find an approximately real solution for each ambiguous locus.

Figure 3 .
Figure 3. Examples of reconstructions for varying noise levels, for a chromosome pair with 60 loci, out of which 50% are ambiguous.Subfigures (a)-(c) show chromosomes simulated with Brownian motion (projected onto the xy-plane), whereas figure (d)-(e) show helix-shaped chromosomes.

Figure 4 .
Figure 4. Comparison between our reconstruction method, ASHIC and PASTIS.The values are the average over 20 runs, with the error bars showing the standard deviation.All experiments took place with 60 loci, with varying levels of noise, as well as varying numbers of ambiguous loci, uniformly randomly distributed over the chromosomes.

Figure 5 .
Figure 5. (a) Reconstruction from a real dataset using our reconstruction method.A dashed line between two beads is used to indicate that there is one or more beads between them, for which we have not given an estimation (due to low contact counts).(b) Bipartite index for the reconstructed chromosomes.The dashed vertical line indicates the known hinge point at locus 146.

Figure 7 .
Figure 7. (a) Total contact counts sorted in increasing order.(b) Ratios between total contact counts.The peak corresponding to the ratio between the 48th and the 47th smallest count is used as a motivation for excluding the 47 loci with smallest total contact from the analysis.(c) Unambiguity quotients for each of the remaining 296 loci, sorted in increasing order.We consider a locus as ambiguous if this ratio is less than 0.4; otherwise, we consider it as unambiguous.

Figure 8 .Figure 9 .
Figure 8. Logarithmic heat maps for the reassigned contact count matrices obtained from the original Patski dataset and from the SNLC reconstruction: (a) and (b) C U ; (c) and (d) C P ; (e) and (f) C A .The axis labels correspond to the 500 unambiguous beads, and the 46 ambiguous loci.